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We use lattice QCD simulations, with MILC gluon configurations and HISQ c-quark propagators, 
to make very precise determinations of moments of charm-quark pseudoscalar, vector and axial- 
vector correlators. These moments are combined with new four-loop results from continuum pertur- 
bation theory to obtain several new determinations of the MS mass of the charm quark and of the MS 
coupling. We find mc(3GeV) = 0.986 (10) GeV, or, equivalently, mdmc) = 1.268 (9) GeV, both for 
n/ = 4 flavors; and aMs(3 GeV, n/ =4) = 0.251 (6), or, equivalently, a^(Mz, n/ = 5) = 0.1174 (12). 
The new mass agrees well with results from continuum analyses of the vector correlator using ex- 
perimental data for e"'"e~ annihilation (instead of using lattice QCD simulations). These lattice and 
continuum results are the most accurate determinations to date of this mass. Ours is also one of 
the most accurate determinations of the QCD coupling by any method. 

PACS numbers: 11.15.Ha,12.38.Aw,12.38.Gc 



INTRODUCTION 



Precise values for the QCD coupling and the 

charm quark's mass rric are important for high-precision 
tests of the Standard Model. Some of the most ac- 
curate mass determinations currently come from zero- 
momentum moments of current-current correlators built 
from the c quark's electromagnetic current (see, for ex- 
ample, Low moments are perturbative and have 
long been known through three- loop order [3, 0, 01 • New 
techniques have recently extended these results to much 
higher moments [1, 0] and, in some cases, to four-loop 
order [1, [l^ . These moments can be estimated non- 
perturbatively, using dispersion relations, from experi- 
mental data for the electron-positron annihilation cross 
section, cr(e+e~ — > 7* X). The c quark's mass is ex- 
tracted by comparing the perturbative and experimental 
determinations. 

In this paper we show how to compute such moments 
directly using accurately tuned, highly realistic numeri- 
cal simulations of QCD in the lattice approximation pdl |. 
As we will show, the correlator moments obtained non- 
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perturbatively from such simulations can be used in place 
of data from e'^e~ annihilation to obtain new, percent- 
accurate determinations of the c quark's mass. With 
lattice QCD, it is also possible to replace the electro- 
magnetic current in the correlator by the pseudoscalar 
operator rricipcl^'^c, thereby providing a completely new 
set of mass determinations and an important cross check 
on the entire methodology. 

In fact the pseudoscalar correlators are particularly 
easy to simulate, and there are no renormalization fac- 
tors required since the corresponding axial- vector current 
is partially conserved in our lattice formalism. Conse- 
quently these correlators give our most accurate masses. 
They also give very accurate values for the QCD cou- 
pling, when combined with our new four-loop results 
from perturbation theory. 

In Section |TT] we describe how to compute pseudoscalar 
correlators and their moments using lattice QCD. We 
discuss techniques for reducing lattice artifacts in Sec- 
tion mil and present new determinations of the c-quark 
mass and QCD coupling from our lattice "data" in Sec- 
tion |TVl In Section El we extend our analysis to include 
vector and axial-vector correlators. We summarize our 
main results in Section IVII In the Appendix we review 
the continuum perturbation theory needed for this anal- 
ysis, including new four-loop results for the pseudoscalar 
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and vector cases. 



II. LATTICE QCD AND PSEUDOSCALAR 
CORRELATORS 



Few-percent accurate QCD simulations have only be- 
come possible in the last few years (see, for exam- 
ple, [13, and accurate simulations of relativistic 
c quarks only in the past year — with the new Highly 
Improved Staggered Quark (HISQ) discretization of the 
quark action [3, [3, which we use here. A lattice QCD 
simulation proceeds in two steps. First the QCD pa- 
rameters — the bare coupling constant and bare quark 
masses in the Lagrangian — must be tuned. Then the 
tuned simulation is used to compute vacuum matrix ele- 
ments of various quantum operators from which physics 
is extracted. An obvious approach to the tuning is to 
choose a lattice spacing a, and then tune each of the 
QCD parameters so that the simulation reproduces the 
experimental value for a corresponding physical quantity 
that is well measured. It is more efficient, however, to 
first choose a value for the bare coupling and then ad- 
just the lattice spacing and bare masses to give physical 
results. 

In the simulations used here, we set the lattice spacing 
to reproduce the correct T' — T meson mass difference 
in the simulations 16], while we tune the u/d, s, c and b 
masses to give correct values for mj, 2rn^— m^, m,,^, and 
mx, respectively. (For efficiency we set m„ = m^; this 
leads to negligible errors in the analysis presented here.) 
The important parameters for the particular simulations 
used in this pap er are listed in TableHJ further details can 
be found in [12, [3 ■ Once these parameters are set, there 
are no further physics parameters, and the simulation will 
accurately reproduce QCD physics for momenta much 
smaller than the ultraviolet (UV) cutoff (A ^ vr/a). 

We have tested these simulations extensively (see, for 
example, [H, [H, [3, [H, Ell ) and, in particular, we have 
done very precise tests for the charm-quark physics most 
relevant to this work. These demonstrate, for example, 
that our simulations reproduce the low-lying spectrum, 
including spin structure, of both charmonium and heavy- 
light mesons {D and Ds) to within our simulation uncer- 
tainties (a few percent or less) Ell ■ 

Given a tuned simulation, it is straightforward to cal- 
culate correlators of the sort used to determine rric- The 
simplest of these is for the c quark's pseudoscalar density, 

G{t)=a' ^(amoe)2(0lj5(x,i)j5(0,0)|0) (1) 



where mpc is the c quark's bare mass (in the lattice La- 
grangian). Here time t is euclidean, and the sum over 
spatial position x sets the total three momentum to zero. 
Note that Git) = G(T - t) = G{T + t) where T is the 
temporal length of the lattice. 



TABLE I: Parameters for the QCD simulations used in this 
paper. The inverse lattice spacing is in units of ri — 
0.321(5) fm [T^ . defined in terms of the static-quark poten- 
tial [13]. L and T are the spatial and temporal size of the 
lattices used for each set of gluon configurations. The con- 
figurations used here were generated by the MILC collabo- 
ration [ij with tt, d and s sea-quarks. The u and d masses 
are set equal to m^/^. The sea-quark masses are given in 
the standard MILC notation which includes a factor of the 
(plaquette) tadpole factor ito. 



Set 


Ti/a 




auQ-mos 


amoc 


Uo 


L/a 


T/a 


1 


2.133(14) 


0.0097 


0.048 


0.850 


0.860 


16 


48 


2 


2.129(12) 


0.0194 


0.048 


0.850 


0.861 


16 


48 


3 


2.632(13) 


0.0050 


0.050 


0.650 


0.868 


24 


64 


4 


2.610(12) 


0.0100 


0.050 


0.660 


0.868 


20 


64 


5 


2.650(8) 


0.0200 


0.050 


0.648 


0.869 


20 


64 


6 


3.684(12) 


0.0062 


0.031 


0.430 


0.878 


28 


96 


7 


3.711(13) 


0.0124 


0.031 


0.427 


0.879 


28 


96 


8 


5.277(16) 


0.0036 


0.018 


0.280 


0.888 


48 


144 



We include two factors of amoc in the definition of 
G{t) so that G{t) becomes independent of the UV cutoff 
as a ^ [isj . Consequently the lattice and continuum 
G{t)s become equal in this limit. Moments G„ are triv- 
ially computed: 



G„ ^ ^(t/a)"G(t), 



(2) 



where, on our periodic lattice [T^, 

t/ae {0,1,2 ... T/2a-l,0, -r/2a + l ...-2,-1}. (3) 
The cutoff independence of G{t) implies that 



Gn 



{amc{fi)y 



+ 0((am,)") (4) 



for n > 4, where md^J.) is the MS mass at scale /i and 
gn is dimensionless. The c mass can be determined from 
moments with n > 6 given G„ from lattice simulations 
and gn from perturbation theory (see Appendix), while 
the QCD coupling can be determined from the dimen- 
sionless moment G4. This assumes that perturbation 
theory is applicable, which should be the case for small 
enough n. 

Note that here and elsewhere in this paper we omit 
annihilation contributions from cc — > gluons — > cc. This 
is allowed provided we omit the same contributions from 
perturbation theory. Annihilation contributions to the 
nonperturbative part of our analysis would be negligi- 
ble (for example, they shift the rj^ mass by approxi- 
mately 2.4MeV, which is less than 0.1% [li.[20|). 
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III. QCD SIMULATIONS 
A. Reduced Moments 

The biggest challenge when using lattice QCD to 
produce c-quark correlator moments is controlling: 1) 
©((aTOc)") errors caused by the lattice approximation; 
and 2) tuning errors in the QCD parameters, and espe- 
cially in the lattice spacing and the c-quark's bare mass. 
We reduce each of these sources of error by making two 
modifications to the moments. 

First we replace G„ by 



Gn 



where Gn^ is the n**^ moment of the correlator to lowest 
order in lattice QCD perturbation theory 21 1, and gi^"* 
is the lowest-order part of gn in continuum perturbation 
theory. The lowest-order on-shell or "pole" mass of the 
c-quark sets the mass scale in the lowest-order lattice 
moments: 
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FIG. 1: Reduced moments i?„ and a ratio of these moments 
from lattice simulations with different lattice spacings a. The 
tight clusters of points at each of the three largest lattice 
spacings correspond to results for different sea-quark masses. 
The dashed lines show the functions used to fit the lattice 
results, with the sea-quark masses set equal to the masses used 
at the smallest lattice spacing. These extrapolation functions 
were used to obtain the a — 0, m„/tj/s — results shown in 
the plot. 



,(0) 



+ 0((amc)") (6) 



In the HISQ formalism, this mass is related to the mass 
moc that appears in the action by [14|: 



.(0) - 



polc.c 



= moc 1 - 



3 (amoc)'' 23 (amoc)^ 



1783 (amoc)' 
537600 



80 2240 
76943 (amoc)^" 
23654400 



(7) 



Introducing Gn removes the explicit factors of the lat- 
tice spacing in the denominator of Eq. ^ , and also can- 
cels finite-a errors to all orders in a and zeroth order 
in as . Thus we expect finite-a errors that are reduced by 
a factor of order as{l/a) w 1/3 when we divide G„ by 
the corresponding lowest-order lattice moment; and we 
find in practice that they are 3-4 times smaller. 

A second modification is to replace the pole mass in 
Gn/Gn'' by the value of the rjc mass obtained from the 
simulation, am^^ (in lattice units) [23 |: 



Gn 



9n 



2am 



(0) 

polc,c / 



(8) 



up to 0{{amc)™'as) corrections. With this additional 
factor, the leading dependence on rricifJ,) enters through 
the ratio mc{fi) / rrij-i^ . Consequently small errors in the 
simulation parameter arriQc are mostly cancelled in this 
expression by corresponding shifts in the simulation value 
for am^^ . This cancellation is accurate up to binding cor- 
rections of order (uc/c)^ w 1/3 in to,,^, and therefore the 



impact of any tuning error in moc is three times smaller 
with this modification [23j . 

Combining these two modifications, we replace G„ by 
a reduced moment: 



G4/Gi°) 



2am 



(0) 

polc.c 



for n = 4, 
Gn/Gl?]'^'^-'^ for.>6, 



The reduced moments can again be written in terms of 
continuum quantities: 

{'^4(aMS' A^/"^c) for n = 4, 
rniaus' m/^c) fQ^^>g^ (10) 
2mc{n)/m^^ ~ 

up to 0{{amc)"^cts) corrections, where r„ is obtained 
from g„ (Eq. ^) and its value, gn \ in lowest-order con- 
tinuum perturbation theory: 




l/(n-4) 



for n — 4, 
for n > 6. 



(11) 



The c mass is obtained from Eq. (jlOp with n > 6 using 
the nonperturbative lattice QCD (LQCD) value for i?„, 
the perturbative QCD (PQCD) estimate for r„, and the 
experimental value for m,,^, 2.980 GeV: 



mdfi) 



7„oxp PQCD 
Vc ' n 



2 pLQCD • 



(12) 



Reduced moment R4 is dimensionless and so depends 
only weakly on mc- Simulation values for this moment 
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can be compared with perturbation theory to obtain esti- 
mates for the QCD couphng: given the c-quark mass and 
a nonperturbative lattice QCD value for R4, we solve the 
equation 

= riiaj^,fi/mc) (13) 

for aj^{fi). Ratios of reduced moments, like Rn/Rn+2 
for n > 6, can also be used in this way to estimate the 
coupling. 



B. Simulation Results 

Our simulation results for Rn{a,mi^/^,ms) are listed 
for different moments n, lattice spacings a, and sea-quark 
masses in Table HIl we also list ratios of reduced moments, 
Rn/Rn+2 for n > 6. Some of these results are plotted ver- 
sus the lattice spacing in Figure [1] Our simulations did 
not include c-quark vacuum polarization, but the correc- 
tion to the moments can be computed using perturbation 
theory (since c-quarks are relatively heavy) |2J] . We find 
that these corrections add 0.7% to R4, and are of order 
0.1% or less for the higher moments considered here. The 
RnS in the table are corrected to include this effect. 

The uncertainty quoted for each Rn{a, m^/^, rris) with 
n > 6 is dominated by the uncertainty in our tuning 
of moc (other sources, such as statistical or finite volume 
errors [2^, are negligible). We tune the bare c-quark 
mass so that our simulations give correct masses for the 
i]c [26| . Our tuning is limited by the precision with which 
we know the lattice spacing a for any given parameter set 
in Table IH since simulations give masses in lattice units 
(that is, arriri^). We determine lattice spacings by com- 
bining MILC's values for ri/a (see Table which are 
accurate to around 0.6%, with a value for ri determined 
from the upsilon spectrum: ri = 0.321(5) fm [l^, which 
is accurate to 1.5%. The corresponding uncertainties in 
the RnS are three times smaller, because of the reduced 
sensitivity of mc/rrijj^ (see above). We include the uncer- 
tainty due to ri/a (0.6%/3 = 0.2%) in the uncertainties 
reported for each separate i?„(a) in Table UIl The un- 
certainty due to ri (1.5%/3 = 0.5%) is included only in 
our final results, after extrapolation (since changes in ri 
affect all i?„(a)s by the same amount). 

Reduced moment i?4 and our ratios of reduced mo- 
ments are much less sensitive to errors in ttiqc because 
the c-mass enters only through radiative corrections, in 
perturbation theory. Consequently uncertainties due to 
mistunings of moc are smaller by an order of magnitude 
or more for these quantities. We account for potential 
tuning errors by assigning an uncertainty of 0.05% to 
each of the R4S and ratios. 

We also include in Table |TT] our results extrapolated to 
zero lattice spacing and zero sea-quark mass p7j . Our 
extrapolation procedure is described in the next section. 



C. arric, mq/mc Extrapolations 

Table HIl and Figure [1] show that our reduced moments 
depend only weakly on the lattice spacing, with most mo- 
ments changing by 0.5% or less between our two smallest 
lattice spacings and only a few percent over our entire 
range. The dependence on the sea-quark masses is even 
weaker. We nevertheless correct our results by fitting the 
variation over our different sets of lattice parameters (Ta- 
ble [H and extrapolating to zero lattice spacing and zero 
sea-quark mass. We do this with a constrained fit '28i. i29j 
of our simulation data to a function of the form: 

Rn{a) =Rn{0) (1 + c„,2(amc)^as + CnA{amc)^as 

+ Cns{amc)*^as + Cn,s{amcfas H ) 

X (l + /„4(2m„/d + ms)/TOcH ) (14) 

where we take mc — 1 GeV and Us — as{\/a). This 
form is motivated by the pattern of a" errors in our lat- 
tice actions, and by chiral perturbation theory, which 
implies nonperturbative corrections that depend linearly 
on the sea-quark masses. (There is sea-quark mass de- 
pendence in perturbation theory, as well, but it enters at 
0{a1{mq/mc)'^) and is neghgible here.) 

The extrapolated results are largely independent of the 
exact functional form used for the extrapolation provided 
reasonable Bayesian priors are included (in the func- 
tion that is minimized in the fit) for each of the coeffi- 
cients Cn^i and fn,i psi. [2^. We use the same Gaussian 
prior, centered at zero with width Uc = 1, for every c„_i, 
for all moments and moment ratios except R4. Moment 
i?4 has larger errors and needs a wider prior; we take 
CTc = 5. We use a prior with width ct/ = 0.1 for the fn,i, 
which is twice as large as the largest coefficient obtained 
from the fits. 

The error estimates from our fits initially increase, and 
decreases, as we add successively higher-order terms 
in the (amc)^ and mq/mc expansions. Eventually the er- 
rors stop increasing and stops changing, again assum- 
ing proper priors for all fit parameters. It is important 
to add terms through this point in order to avoid under- 
estimating uncertainties in the final fit results. Adding 
further terms has no effect on fit results (values or errors). 

Only a single term is needed in each series to get good 
fits for Rn with n > 6 if we discard data from the largest 
lattice spacing; and our final results (values and errors) 
are little changed. We, however, retain results from the 
coarsest lattices, despite the large value of amc for those 
lattices, in order to test our priors. Fitting all of our 
simulation data, we get good fits with two terms in the 
(aTOc)^ and a single term in the mq/mc expansion. To be 
certain of convergence, we used eight terms in the first 
expansion and two in the second to obtain the extrapo- 
lated results in Table ITIl The fact that we get good fits 
(x^ per data point less than one) even when we include 
data from the coarsest lattices helps validate the design 
of our fit function and priors, and it reassures us that our 
fits are not underestimating errors. 
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TABLE II: Simulation results for Rn(a,m^/^,nis) for different lattice parameter sets (see Table |T]|. The inverse lattice spac- 
ing is in GeV. Extrapolations to zero lattice spacing and zero sea-quark masses are given for each quantity, together with 
the corresponding value for mc(/i) (in GeV) or aMs(pi) for n/ — A flavors and = 3 GeV. 



Set: 1 


2 


3 


4 


5 


6 


7 


8 






a-^: 1.31 


1.31 


1.62 


1.60 


1.63 


2.26 


2.28 


3.24 






Re, 1.448(3) 


1.447(3) 


1.494(3) 


1.492(3) 


1.491(3) 


1.514(3) 


1.511(3) 


1.519(3) 


1.528(11) 


0.986(10) 


Rs 1.372(3) 


1.371(3) 


1.387(3) 


1.386(3) 


1.384(3) 


1.374(3) 


1.373(3) 


1.370(3) 


1.370(10) 


0.986(11) 


Rio 1.329(3) 


1.328(3) 


1.326(3) 


1.326(3) 


1.324(3) 


1.306(3) 


1.305(3) 


1.304(3) 


1.304(9) 


0.973(19) 


Ri2 1.294(3) 


1.293(3) 


1.284(3) 


1.284(3) 


1.281(3) 


1.263(3) 


1.262(3) 


1.262(3) 


1.265(9) 


0.969(23) 


Ri4 1.264(3) 


1.264(3) 


1.252(2) 


1.251(2) 


1.248(2) 


1.232(2) 


1.231(2) 


1.232(2) 


1.237(9) 


0.967(28) 


Ri6 1.239(2) 


1.239(2) 


1.228(2) 


1.226(2) 


1.223(2) 


1.207(2) 


1.206(2) 


1.210(2) 


1.215(9) 


0.965(33) 


Rw 1.218(2) 


1.218(2) 


1.208(2) 


1.205(2) 


1.202(2) 


1.187(2) 


1.187(2) 


1.191(2) 


1.198(9) 


0.963(38) 




Set: 1 


2 


3 


4 


5 


6 


7 


8 






a-^: 1.31 


1.31 


1.62 


1.60 


1.63 


2.26 


2.28 


3.24 




"ms(m) 


i?4 1.162(1 


) 1.161(1) 


1.189(1) 


1.187(1) 


1.187(1) 


1.223(1) 


1.221(1) 


1.249(1) 


1.281(5) 


0.252(6) 


Rq/Rs, 1.055(1 


) 1.055(1) 


1.078(1) 


1.076(1) 


1.077(1) 


1.101(1) 


1.101(1) 


1.109(1) 


1.113(2) 


0.249(6) 


Rs./Rw 1.033(1 


) 1.033(1) 


1.046(1) 


1.045(1) 


1.046(1) 


1.052(1) 


1.052(1) 


1.051(1) 


1.049(2) 


0.224(31) 


R10/R12 1.027(1 


) 1.027(1) 


1.033(1) 


1.033(1) 


1.034(1) 


1.034(1) 


1.034(1) 


1.033(1) 


1.031(2) 


0.241(30) 


R12/R14, 1.023(1 


) 1.023(1) 


1.025(1) 


1.026(1) 


1.026(1) 


1.025(1) 


1.025(1) 


1.024(1) 


1.022(2) 


0.243(47) 


Ri4/Rie 1.020(1 


) 1.020(1) 


1.020(1) 


1.021(1) 


1.021(1) 


1.020(1) 


1.020(1) 


1.019(1) 


1.017(2) 


0.242(70) 


i?16/i?18 1.017(1 


) 1.017(1) 


1.016(1) 


1.017(1) 


1.017(1) 


1.017(1) 


1.017(1) 


1.016(1) 


1.014(2) 


0.241(96) 



Moment R4 and the ratios of moments are more accu- 
rately determined in our simulation than the other i?„s, 
and so typically require an additional term in the {arric)^ 
expansion. Again, however, the eight terms we use are 
many more than the minimum needed. 

Our final error estimates depend upon the widths of 
our priors [1^. We tested these widths in a couple of 
ways, beyond including simulation data from the coars- 
est lattices. First we compared our widths with the val- 
ues suggested by the empirical Bayes procedure described 
in [2^. This procedure uses the variation in the data it- 
self to determine, for example, an optimal value for CTc. 
The widths we use are two to four times larger that what 
is indicated by the empirical Bayes criterion, suggesting 
that our error estimates are conservative. The dominant 
fit coefficients in the (awc)^ expansion for /Jgi for ex- 
ample, range between —0.05 and —0.20, which is much 
smaller than the CTc = 1 we use. 

As a second test, we verified that our extrapolation 
procedure gives consistent results when data from either 
the smallest or the largest lattice spacing is discarded. 
That is, we demonstrated that results obtained from the 
truncated data sets agree within errors with results from 
the full set of simulation data. This shows that our error 
estimates are robust even when working with limited sim- 
ulation data sets. As mentioned above, our final results 
are not much affected by data from the coarsest lattice 
spacing. Simulation data from the finest lattice spacing, 
on the other hand, has a very significant impact. 
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FIG. 2: mc{fi), for fi = 3 GeV and n/ = 4 flavors, from dif- 
ferent moments of correlators built from four different lattice 
operators. The gray band is our final result for the mass, 
0.986 (10) GeV, which comes from the first two moments of 
the pseudoscalar correlator (upper- left panel). 



IV. EXTRACTING mc(fi) AND a^{fi) 



To convert the extrapolated reduced moments into 
c masses and coupling constants, we require perturba- 
tive expansions for the r„ in Eq. (fT^ . These are easil; 
computed from the expansions for gn [1, 0, [E B S B 
using Eq. (111^ : details can be found in the Appendix. 
The perturbative expansions have the form 



rn = l+r„^laJj^{^l)+rn,2a^{^J.)+rn,3aJ^{^l) + . . . (15) 
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FIG. 3: ajj^g{Mz;nf — 5) from R4 and ratios Rn/Rn+2- The 
gray band is our final result for the coupling, 0.1174(12), 
which comes from R4 and Ra/Rs- 



where we set the renormalization scale /i to 3GeV [30|. 
The full third-order coefficients for the n — 4,6,8 mo- 
ments were computed for this analysis and are presented 
in the Appendix. The third-order coefficients for mo- 
ments with n > 10 are only partially complete: our anal- 
ysis includes all /x-dependent terms (that is, log"(/i/TOc) 
terms), but the constant parts have not yet been com- 
puted. Consequently we take the truncation uncertainty 
in r„ to be of order [sH 

fC^«|ls('^) forn = 4,6,8, 
<^r^ = <^ (16) 
[r^-al^if,) forn>10, 

where 

C- = max(|r„,i|,|r„,2Ur„,3|). (17) 

Another source of uncertainty in all of our moments 
comes from nonperturbative effects. In the previous sec- 
tion, we discuss how we remove nonperturbative con- 
tributions involving the sea-quark masses. To assess 
the importance of gluonic contributions, we also include 
the leadi ng g luon-condensate contribution in our mo- 
ments 0, IS^TISSj. We do this by multiplying by a 
factor of the form {l + dn{asG^ /tt) /{2mc)*)) where, here, 
rric — mc{mc) and d„ is computed through leading order 
in ay^{mc). The value of the condensate is not well 
known; we set {asG^/ir) = ± 0.012 GeV", which covers 
the range of most current estimates [34]. 

Note that coefficients in the r„ expansion, Eq. (fT5|l . 
depend upon mj,(/i) through scale-dependent loga- 
rithms, log"(/i/r7ic(/i)). Consequently, the mass appears 
on both sides of Eq. (fT2)) . and the equation is an im- 
plicit equation for mc{^J■)■ The rn^ (M)"dependence on the 
right-hand side, however, is suppressed by a-^di), and 
therefore the equation is easily solved numerically. 

Our final results for the c-quark's mass, rndp) at 
/i = 3GeV for Uf = A flavors in the MS scheme, are 
listed in Table |TT1 and plotted in the upper-left panel of 
Figure [D As is clear from the figure, all moments agree 



TABLE III: Sources of uncertainty in the determinations of 
mc{i-L = 3GeV,n/ =4) and ajjg{Mz,nf — 5) from different 
reduced moments i?„ of the pseudoscalar correlator. The un- 
certainties listed are percentages of the final results. 





mc(/i) 




Mb 


(Mz) 




Re 


Ra 


R4 


Rg/Rs 


extrapolation 


0.2% 


0.3% 


0.4% 


0.2% 


perturbation theory 


0.4 


0.3 


0.6 


0.6 


a^jjg- uncertainty 


0.3 


0.4 






'mc{fJ.) uncertainty 






0.1 


0.1 


gluon condensate 


0.3 


0.0 


0.4 


0.7 


statistical errors 


0.1 


0.0 


0.2 


0.1 


moc errors from ri/a 


0.5 


0.6 


0.3 


0.4 


moc errors from ri 


0.6 


0.6 


0.1 


0.1 


^u/d/s extrapolation 


0.2 


0.2 


0.2 


0.4 


finite volume 


0.1 


0.1 


0.0 


0.3 


fj, — > Mz evolution 


0.0 


0.0 


0.1 


0.1 


Total 


1.0% 


1.1% 


1.0% 


1.1% 



on the mass although the higher moments may be less 
trustworthy (see 0). The first two moments {n = 6,8) 
give results that are twice as accurate as the others be- 
cause we have full 0(ajj^) perturbation theory in these 
cases. We average the two results, which agree, to obtain 
our final result for the mass: 

mc(3GeV,n/=4) = 0.986(10)GeV (18) 

Evolving down to scale /i = mc{fi) using fourth-order 
evolution [H, [H, [13, HI] , this is equivalent to [SOj 

TOc(mc, n/ =4) = 1.268 (9) GeV. (19) 

We used a^jg (3 GeV, n/ = 4) = 0.252 (10) in the per- 
turbation theory needed to extract mc{fi). We derived 
this from the current Particle Data Group average for 
the 71/ = 5 coupling at ^ = M^, which is 0.1176 (20) 
The coupling can also be extracted directly from R4 and 
from the ratios Rn/Rn+2, as discussed above. Taking 
rricip) — 0.986 (10) GeV, we obtain the couplings, for 
scale ^ = 3 GeV and Uf — 4, shown in Table [ill The 
first two determinations listed in the table are far more 
accurate than the others because we know perturbation 
theory through third order. We can average these to ob- 
tain a composite value for the coupling of: 

aj;js(3GeV,n/ = 4) = 0.251 (6). (20) 

To allow comparison with other work we converted 
our couplings to = 5 by adding a 6-quark with 
mass mb{mb) = 4.20(7) GeV [40|, and evolving them to 
scale Mz ■ The results are shown in Figure [3l Averaging 
the first two numbers, which agree with each other, we 
get: 

aMs(Mz,n/ = 5) = 0.1174(12). (21) 



7 



The leading sources of uncertainty in mc{n) and 
ay^{Mz) are listed in Table IIIII for those calcula- 
tions where we have full perturbation theory through 
0{a^j^) [2^. The dominant uncertainty in the masses 
comes from potential tuning errors in the c-quark masses 
used in the simulation. Truncation errors from pertur- 
bation theory dominate for the coupling, with nonper- 
turbative contributions from the gluon condensate also 
becoming important. In addition to the various sources 
discussed above, there are also uncertainties due to the 
finite spatial volume of our lattices; our lattices were ap- 
proximately 2.5 fm across. While our simulations showed 
no measurable volume dependence [1^ , lattice perturba- 
tion theory shows finite- volume sensitivity for the higher 
(more infrared) moments. This is negligible for lower mo- 
ments but grows with n. The finite-volume sensitivity is 
mostly an artifact of perturbation theory; confinement 
significantly reduces finite- volume effects. Consequently 
we assign a finite-volume error to our perturbative fac- 
tors that is equal to the entire finite- volume correction in 
perturbation theory. 



V. mdp) FROM OTHER CORRELATORS 

The close agreement on rric between different moments 
is important evidence that we understand our systematic 
errors since these enter quite differently in different mo- 
ments. To further check this we repeated our analysis 
for three different correlators, which we formed by re- 
placing the pseudoscalar operator mocjs with each of the 
following c-quark currents on the lattice: 



(1) 



Jr = V'c(a;)7MV'c(x), 



(22) 
(23) 
(24) 



The first two currents are different lattice discretizations 
of the vector current and were evaluated for space- like fis; 
and the first of these was evaluated in Coulomb gauge. 
The third current is a lattice discretization of the axial 
vector current and was evaluated for time- like fi. The 
superscript on each j labels the "taste" carried by that 
operator, using the notation presented in the Appendices 
of Taste is a spurious quantum number, analogous 
to flavor, that is an artifact of staggered-quark lattice 
discretizations like the HISQ formalism. Taste should not 
affect physical results and therefore operators carrying 
different taste here should give identical results in the 
a — > limit. By studying these different currents, we 
not only test for conventional systematic errors, but also 
verify that HISQ-specific taste effects are negligible [4lj . 

A complication in our lattice analysis of these vector 
(or axial- vector) correlators is that none of the currents is 
conserved (or partially conserved) on the lattice. Conse- 
quently, each lattice current is related to its correspond- 



TABLE IV; Simulation results for the reduced moments Rn\ 
extrapolated to a = 0, from correlators of local axial-vector 
and vector lattice currents, and a point-split lattice vector 
current. Corresponding values for mc(/x) (in GeV), for fi = 
3 GeV and n/ = 4, are also given. Only results for parameter 
sets 1, 4 and 6 from Table |I] were used for the first and last 
currents; results from these sets were combined with results 
from set 8 (the smallest lattice spacing) for . 



,-(5a.) 



R 



R 



(3) 



,•(1) 



10 
12 
14 
16 
18 



1.240(27) 0.97(3) 

1.159(25) 0.97(5) 

1.126(24) 0.99(5) 

1.103(24) 0.99(6) 

1.082(23) 0.99(8) 

1.064(23) 1.00(9) 



1.233(16) 1.00(3) 

1.183(15) 1.00(4) 

1.162(15) 0.98(5) 

1.139(15) 0.97(6) 

1.120(15) 0.97(8) 

1.106(14) 0.97(10) 

1.093(14) 0.97(13) 



1.261(30) 
1.163(27) 
1.132(27) 
1.113(26) 
1.094(26) 
1.076(25) 
1.059(25) 



0.97(4) 
1.02(5) 
1.02(6) 
1.01(7) 
1.01(9) 
1.01(11) 
1.02(14) 



ig continuum operator by a renormalization constant: 



Zij) = ^0)(Q,_(7r/a),amoc) 
where j is one of the lattice currents \ 



(25) 



and jcont is the continuum current = V'T/jV' for the 
first two js and = V'757mV' for the last. Consequently 
moments of the correlators of these lattice currents have 
the form 



1 gi^''°"'\aj^{ij),ii/mc) 



(26) 



where Z^^^'^Gn^ is the continuum result for n > A. To 
cancel the renormalization factor we redefine the reduced 
moments for these correlators to be 



7O) 



1/2 



-.00) 



„Ocont) 



2?nc(/i)/mO) 



(27) 
(28) 



where n > 6, and to'-') is the tp mass for the vector cur- 
rents (which couple to the V') and the 7]^ mass for the 
axial- vector current. Again we divide each moment Gn •* 
by its value g1''°^ in leading-order lattice perturbation 
theory in order to minimize finite-lattice-spacing errors. 
And again the perturbative expansion for ri"'°°°*'' can be 
obtained from continuum perturbation theory expansions 

for the gn°°'^^^ (see the Appendix). 

(7) 

Our simulation results for i?„ , extrapolated to lat- 
tice spacing a = 0, are given in Table IIVI for different 
moments n and each of the three currents [4^ . Per- 
turbative coefficients for the vector-current rn''°"*^s are 
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1.05 




0.95 



10 



n 



FIG. 4: Ratio of the extrapolated simulation results from 



Table HVI for with j 



to results derived from ex- 



periment in [3 for different moments n. 



discussed in the Appendix; the coefficients for the tem- 
poral axial-vector current can be derived from the pseu- 
doscalar coefficients (also in the Appendix) using Ward 
identities 

By combining perturbative with nonperturbative re- 
sults, we obtain the values for mdn), with /i = 3GeV 
and n/ = 4, that are listed in Table [TVl and plotted in the 
top-right and bottom panels of Figure [H Results from 
all moments agree with each other and with the pseu- 
doscalar result (the gray band in the plots), although 
here the errors are about twice as large for the smaller 
moments. 

Values for the lowest four moments of the vector corre- 
lator are derived from experimental data for e^e~ annihi- 
lation in Converting these into reduced moments, we 

compare them with our extrapolated Rn^s for j = jl^^ 
(from Table lIVp in Figure HI We find that experiment 
and our simulation results agree to within combined er- 
rors of better than 2%. This comparison is more accurate 
than comparing masses, because we are comparing the 
(reduced) moments directly, without recourse to further 
perturbation theory. 



VI. CONCLUSIONS 

Our results are by far the most accurate determination 
of the c-quark mass from lattice QCD [31 . Such precision 
is possible because the matching between lattice param- 
eters and continuum parameters here relies upon contin- 
uum perturbation theory, which is much simpler than 
lattice QCD perturbation theory. Consequently pertur- 
bation theory can be pushed to much higher orders. The 
precision of this continuum calculation is matched by 
that of our nonperturbative lattice analysis because of 
the quality of the MILC configurations and of our highly 
corrected HISQ action for c-quarks. 

The agreement between masses from different mo- 
ments, and from different correlators — 27 determina- 
tions in all — is an important check on systematic errors 
of all sorts since these enter in very different ways in each 



> 1-2 

a; 

O 1.1 

3 1 

u 

S 0.9 

> 1-2 
C5 1.1 

3 1 

S 0.9 



.(5) 
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51 



20 24 28 32 36 40 44 48 52 56 60 
Moment n 

FIG. 5: rricip. — 3GeV) from large-n moments of the pseu- 
doscalar and the (local) vector correlators. The gray band is 
our final result for the mass (from n = 6, 8), 0.986 (10) GeV. 
The perturbative part of the analysis was evaluated at fi = 
nicip) using formulas from Q, and the results evolved to 
= 3 GeV using fourth-order evolution. Uncertainties due 
to the gluon condensate are not included (see text). 



calculation. Note that the different reduced moments in 
our analysis vary in value by as much as 43% (from 1.06 
to 1.52), and yet they all agree on the value of mc to 
within a few percent once we account for differences in 
the perturbative parts. 

One surprising feature of our results is that even the 
higher moments give correct values for the quark mass, 
albeit with larger errors. Nonperturbative effects grow 
with n but our results show no systematic deviation un- 
til very large n, as is evident in Fig. O which shows results 
for 20 < n < 62. We have not included potential errors 
due to the gluon condensate in this figure. The error 
bars would have been much larger had we done so. For 
example, they would have been about 5 times larger at 
n = 40 in the pseudoscalar plot (16% rather than 3%). 
This might suggest that the condensate is smaller than 
we allowed for — say {a.G^/n) < 0.003 GeV — but we 
have not analyzed this carefully enough to make a strong 
statement. The error bars shown in the plots start to 
grow rapidly just where it becomes clear that perturba- 
tion theory is failing (because of large coefficients). 

Our lattice result for the mass, TOc(3GeV,n/ = 
4) = 0.986 (10) GeV, agrees well with the continuum de- 
termination from e'^e" annihilation data, which gives 
0.986 (13) GeV This provides further strong evidence 
that the different systematic errors in each calculation 
are understood. Similarly our new value for the cou- 
pling, ay^{AIz,nf — 5) = 0.1174 (12), agrees very well 
with non-lattice determinations kol j45| and our other 
determinations from lattice QCD [IJ, [2^ . It is also more 
accurate than most determinations. 

The close agreement of our results with non-lattice de- 
terminations of the mass and coupling, and the taste- 
independence of our masses, is also further evidence that 
the simulation methods we use are valid. While early con- 
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cerns about the light-quark discretization used here have 
been largely addressed [H, \^ , it remains important to 
test the simulation technology of lattice QCD at increas- 
ing levels of precision, given the critical importance of 
lattice results for phenomenology. 

Our results are particularly relevant to the recent, very 
accurate analysis of tt, D and meson decay con- 
stants using the same HISQ formalism for valence quarks 
and many of the same MILC configuration sets that we 



use here [15| . The correlators in our pseudoscalar analysis 
are identical to those used to extract the decay constants 
in the earlier study, except that here the light valence 
quarks have been replaced by c quarks, which should 
make finite-a errors worse. The agreement of our pseu- 
doscalar analysis of Rg with the continuum analysis for 
rricifJ-) indicates that our extrapolated lattice cc correla- 
tors are reliable to within 2% or better. The accuracy of 
our cc results, together with the demonstrated accuracy 
of the TT and K results in [T5| from light-quark corre- 
lators, strongly suggest that the corresponding D and 
Ds predictions, from correlators with a light quark and 
a c quark, are reliable. This conclusion is made more 
important by the 3.6f7 discrepancy between the pre- 
diction and recent experimental results jH, |4^, [s^l . 

The lattice analysis will be improved as data becomes 
available for smaller lattice spacings. Also a very accu- 



rate c-quark mass will allow us to make similarly accurate 
determinations of the s-quark mass 5l[ . This is because 
the ratio ms/rric can be determined very accurately in 
lattice simulations where the s and c quarks are analyzed 
using the same formalism, as here. Finally four-loop per- 
turbation theory for additional moments would improve 
both lattice and continuum determinations. 
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^(33) 


1 


1.3333 


3.1111 


0.0000 


0.1154 


-6.4815 


0.0000 


-1.2224 


2.5008 


13.5031 


0.0000 


2 


0.5333 


2.0642 


1.0667 


7.2362 


1.5909 


-0.0444 


7.0659 


-7.5852 


0.5505 


0.0321 


3 


0.3048 


1.2117 


1.2190 


5.9992 


4.3373 


1.1683 


14.5789 


7.3626 


4.2523 


-0.0649 


4 


0.2032 


0.7128 


1.2190 


4.2670 


4.8064 


2.3873 




14.7645 


11.0345 


1.4589 


5 


0.1478 


0.4013 


1.1821 


2.9149 


4.3282 


3.4971 




16.0798 


16.6772 


4.4685 


6 


0.1137 


0.1944 


1.1366 


1.9656 


3.4173 


4.4992 




14.1098 


19.9049 


8.7485 


7 


0.0909 


0.0500 


1.0912 


1.3353 


2.2995 


5.4104 




10.7755 


20.3500 


14.1272 


8 


0.0749 


-0.0545 


1.0484 


0.9453 


1.0837 


6.2466 




7.2863 


17.9597 


20.4750 


1 


1.0667 


2.5547 


2.1333 


2.4967 


3.3130 


-0.0889 


-5.6404 


4.0669 


0.9590 


0.0642 


2 


0.4571 


1.1096 


1.8286 


2.7770 


5.1489 


1.7524 


-3.4937 


6.7216 


6.4916 


-0.0974 


3 


0.2709 


0.5194 


1.6254 


1.6388 


4.7207 


3.1831 




7.5736 


13.1654 


1.9452 


4 


0.1847 


0.2031 


1.4776 


0.7956 


3.6440 


4.3713 




4.9487 


17.4612 


5.5856 


5 


0.1364 


0.0106 


1.3640 


0.2781 


2.3385 


5.3990 




0.9026 


18.7458 


10.4981 


6 


0.1061 


-0.1158 


1.2730 


0.0070 


0.9553 


6.3121 




-3.1990 


16.9759 


16.4817 


7 


0.0856 


-0.2033 


1.1982 


-0.0860 


-0.4423 


7.1390 




-6.5399 


12.2613 


23.4000 


8 


0.0709 


-0.2660 


1.1351 


-0.0496 


-1.8261 


7.8984 




-8.6310 


4.7480 


31.1546 



TABLE V: Moments of the pseudoscalar (upper 8 lines) and the vector (lower 8 lines) correlators, where currently unknown 
coefficients are denoted by a dash and set to zero in our analysis. The numbers correspond to QCD with one massive c-quark 
and three massless {u, d, s) quarks. 



Appendix A: Continuum perturbation theory 



by 



The correlators of two pseudoscalar {itpcjdtpc), vector 
ipcln'4^c^ and axial-vector ipdnlb^c currents are defined 



q^W{q^) = I / da;e^'?^(0|rj5(x)j5(0)|0), (29) 
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{q^qu - q'g^.Wiq') + q^q.Uiiq') = n%{q) 

= 1 jdxe^'^^{)\Tjl{x)jtm{)). (30) 

where S — v and a for the vector and axial-vector cases, 
respectively. The polarization functions can be expanded 
in the variable z = q^ / {2mc{iJi)Y : 



167r2 



(31) 



The polarization function and the c-quark mass mc{ii) 
are renormalizcd in the MS-scheme. 

The longitudinal part of the axial-vector current 
(which is of interest in the present context) and the 
pseudoscalar correlator are related by the axial Ward- 
identity ^Si]: 

qf^q^ni^^iq) = {2m)^q^nP (q^) + contact term. (32) 
Comparing different orders in z, 



where Im^ = log(m^(/x)//x^). The coefficients cj^''^ up 
through fc = 8 are listed in Table fVl for both pseudoscalar 
and vector correlators. 



The four-loop coefficients C'^"'' and C2^^'' for the pseu- 
doscalar correlator are ne wjssj : and c['^°^ for the vector 
correlator comes from [J ^The four-loop coefficients 

^3^°^ for the pseudoscalar case is also new [l^j, as is 

6*2^°^ for the vector case [131 • The order terms up 
through k = 8 are given in while results for 

higher /c- values are given in [a Mi See also jsH^ for nf- 
dependent four-loop terms, and [7| for the pseudoscalar 
case. Throughout this paper we take the number of light 
(massless) active quark flavors to be n; = 3, and the 
number of heavy (massive) quarks is set to n^, = 1. For 
numerical work, we use fi — 3 GeV. 

The expansion coefficients Ck are related to the coef- 
ficients 5„ of Eq. ^ by: 



^(30) 



fc+1 — 



for fc > — 1, 



(33) 



allows us to extract moments of the pseudoscalar correla- 
tor from those of the longitudinal part of the axial- vector 
correlator, and vice versa. The contact term contributes 
only to k = —2. The coefficients of the perturbative ex- 
pansion depend logarithmically on rric and can be written 
in the form 



92k+2 

(0) 
92k+2 



Ck 



C. 



(0) 



(35) 



and the coefficients rn,i in Eq. p5|) are obtained through 
the series expansion of Eq. (fTTj) . For the vector corre- 



lator the coefficients 
expansions of 



„0m) 



are defined through the series 



(0) , asin) /^(lo) , ^(11) 



-I- 



(20) 



^k " 



^(30) ^(31), 
^k ' ^k " 



, /=.(32),2 
+'-^k *m 



/=.(33),3 



.(34) 



' 2fe+2 



cv 



1/2 



(36) 
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